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pg ■ Abstract. This paper reviews and extends some recent results on the multivariate frac- 

. , tional Brownian motion (mfBm) and its increment process. A characterization of the 

f^ f mfBm through its covariance function is obtained. Similarly, the correlation and spectral 

■^ I analyses of the increments are investigated. On the other hand we show that (almost) 

■ all mfBni's may be reached as the limit of partial sums of (super)linear processes. Fi- 

^vj , nally, an algorithm to perfectly simulate the mfBm is presented and illustrated by some 

simulations. 

^'- 1. Introduction 

-)— » ■ 

The fractional Brownian motion is the unique Gaussian self-similar process with station- 
ary increments. In the seminal paper of Mandelbrot and Van Ness |22], many properties 
of the fBni and its increments are developed (see also [30] for a review of the basic proper- 
ties). Depending on the scaling factor (called Hurst parameter), the increment process may 
00 ! exhibit long-range dependence, and is commonly used in modeling physical phenomena. 

^ I However in many fields of applications (e.g. neuroscience, economy, sociology, physics, 

Q ■ etc), multivariate measurements are performed and they involve specific properties such as 

f>^ . fractality, long-range dependence, self-similarity, etc. Examples can be found in economic 

O I time series (see jTT], [M], |l5]), genetic sequences p], multipoint velocity measurements in 

^-^ ■ turbulence, functional Magnetic Resonance Imaging of several regions of the brain [1]. 

It seems therefore natural to extend the fBni to a multivariate framework. Recently, this 
question has been investigated in f2U\ [T21 [3]. The aim of this paper is to summarize and to 
r> ! complete some of these advances on the multivariate fractional Brownian motion and its 

C^ \ increments. A multivariate extension of the fractional Brownian motion can be stated as 

follows : 

Definition 1. A Multivariate fractional Brownian motion (p-mfBm or mfBm) with pa- 
rameter H G (0, ly is a p -multivariate process satisfying the three following properties 

• Gaussianity, 



(N 
> 



1991 Mathematics Subject Classification. 26A16, 28A80, 42C40. 

Key words and phrases. Self similarity ; Multivariate process ; Long-range dependence ; Superlinear 
process ; Increment process ; Limit theorem. 

Research supported in part by ANR STARAC grant and a fellowship from Region Rhone- Alpes (France) 
and a Marie-Curie International Outgoing Fellowship from the European Community. 
Research supported in part by ANR InfoNetComaBrain grant. 

1 



2 AMBLARD, COEURJOLLY, LAVANCIER, AND PHILIPPE 

• Self- similarity with parameter H G (0, 1)^, 

• Stationarity of the increments. 

Here, self-similarity has to be understood as joint self-similarity. More formally, we use 
the following definition. 

Definition 2. A multivariate process {X(t) = (Xi(t), ■ ■ ■ , Xp(t)))t(zK is said self-similar if 
there exists a vector H = {Hi, ■ • • , Hp) G (0, 1)^ such that for any A > 0, 

(1) (Xi(At),--- ,X,(At))),eR = (A^^Xi(t),--- ,A^-X,(t))^^^, 

where = denotes the equality of finite- dimensional distributions. The parameter H is called 
the self- similarity parameter. 

This definition can be viewed as a particular case of operator self-similar processes by 
taking diagonal operators (see [T2l [T6l [T71 |2T]). 

Note that, as in the univariate case [12], the Lamperti transformation induces an isometry 
between the self-similar and the stationary multivariate processes. Indeed, from Definition 
[21 it is not difficult to check that (F (t))igK is a p-multivariate stationary process if and only 
if there exists H e (0, 1)^ such that its Lamperti transformation (t^^Yi(\og{t)) , . . . , t^pl^(log(t)))j(=iR 
is a p-multivariate self-similar process. 

The paper is organized as follows. In Section |2l we study the covariance structure of 
the mfBni and its increments. The cross-covariance and the cross-spectral density of the 
increments lead to interesting long-memory type properties. Section [3] contains the time 
domain as well as the spectral domain stochastic integral representations of the mfBm. 
Thanks to these results we obtain a characterization of the mfBm through its covariance 
matrix function. Section|3]is devoted to limit theorems, the mfBm is obtained as the limit of 
partial sums of linear processes. Finally, we discuss in Section [5] the problem of simulating 
sample paths of the mfBm. We propose to use the Wood and Chan's algorithm |32] 
well adapted to generate multivariate stationary Gaussian random fields with prescribed 
covariance matrix function. 

2. Dependence structure of the mfBm and of its increments 

2.1. Covariance function of the mfBm. In this part, we present the form of the co- 
variance matrix of the mfBm. 

Firstly, as each component is a fractional brownian motion, the covariance function of 
the i-th component is well-known and we have 

2 

(2) EX^{s)Xi{t) = ^ {\s\^^^ + \t\^^^ - \t - s\^"'} . 

with af := var(Xj(l)). The cross covariances are given in the following proposition. 

Proposition 3 (Lavancier et al. [20|). The cross covariances of the mfBm satisfy the 
following representation, for all {i,j) G {1, . . . ,p}^, i 7^ j , 
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[1) IfHi+Hj 7^ 1, there exists {pij,r]ij) G [—1, l]xM with pij = pj^i = corr(Xi(l),Xj(l)) 
and rjij = —rjj^i such that 



(p„-r/,,,sign(t-s))|t-s|^»+^^}. 



(3) EX,(s)X,(t) = ^ {(p,j + r/,,,sign(s))|s|^»+^^ + (p,,, - r/,jSign(t))|ti-»--^ 



(2) IfHi+Hj = 1, there exists {pij,r]ij) G [—1, l]xM with pij = pj^i = corr(Xj(l),Xj(l)) 
and ?7j j = —fjji such that 
(4) 
EXi{s)Xj{t) = ^ {p^j(|s| + |t| -\s- t\) + fii,,it\og \t\ -slog \s\ -(t-s) log \t - s\)} . 

Proof. Under some conditions of regularity, Lavancier et al. [20] actually prove that Propo- 
sition 3 is true for any L^ self-similar multivariate process with stationary increments. The 
form of cross covariances is obtained as the unique solution of a functional equation. Formu- 
lae ([3]) and (j4]) correspond to expressions given in |20] after the following reparameterization 
• Pi,j ~ i^ij + ^j,i)/'^ ^^^ Vij = (cjj- — Cj,j)/2 where Cjj and Cj^i arise in [2U] . 

D 

Remark 1. Extending the definition of parameters pijjpijjrjijjfjij to the case i = j, we 
have pi^i = pi^i = 1 and rji^i = fji^i = 0, so that ([2]) coincides with Q and (j4]). 

Remark 2. The constraints on coefficients pij, pij, rjij, fjij are necessary but not sufficient 
conditions to ensure that the functions defined by ([3]) and (^ are covariance functions. This 



problem will be discussed in Section 3.4 



Remark 3. Note that coefficients pij, pij,T]ij,fjij depend on the parameters {Hi,Hj). As- 
suming the continuity of the cross covariances function with respect to the parameters 
{Hi,Hj), the expression (jl]) can be deduced from ([3]) by taking the limit as Hi + Hj tends 
to 1, noting that ((s + 1)-^ - s^ - 1)/(1 ~ H) ^ slog|s| - (s + l)log|s + 1| as H ^ I. 
We obtain the following relations between the coefficients : as Hi + if j — )■ 1 

Pi,j ~ Pi,j and {I- Hi- Hj)r]ij ~ fjij. 

This convergence result can suggest a reparameterization of coefficients rjij in [1 — Hi — 

2.2. The increments process. This part aims at exploring the covariance structure of 
the increments of size 5 of a multivariate fractional Brownian motion given by Definition [1] 
Let AsX = {X{t+6)—X(t))t^^ denotes the increment process of the multivariate fractional 
Brownian motion of size 6 and let AgXi be its i-th component. 

Let '~fij{h,6) = E,AsXi{t)AsXj{t + h) denotes the cross- covariance of the increments of 
size S of the components i and j. Let us introduce the function WijQi) given by 

_ r {p,, - r/,,,sign(/i))|/i|^'+^^ if H, + Hj ^ 1, 
^ ' '''^ ^ \ hj\h\ + m,h\og \h\ if H, + H^ = I. 
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Then from Proposition |3l we deduce that '-fij{h,6) is given by 



(^iCTj 



(6) 7ij(^, S) = -Y' \ ^».i(^ -5)- 2wij{h) + Wij{h + 6)j. 

Now, let us present the asymptotic behaviour of the cross-covariance function. 
Proposition 4. As \h\ — )■ +oo, we have for any 6 > 

(7) ^i,j{h, 6) ~ a,a,5>|^'+^^-2Ki,,(sign(/i)), 
with 

(8j /.,, (sign(/zj j - I ^^^^.gigj,^^) ,^ ^^ + ^^. ^ i_ 

Proof. Let a = Hi+Hj. Let us choose /?,, such that |/i| > 5, which ensures that sign{h — 6) 
sign(/i) = sign(/i + (5). When a ^ 1, this allows us to write 



^i^Uula 



li,j{K ^) = -i^\K {Pi,j - ^ijsign(/i)) B{h), 

with B{h) = (1 - f)" - 2 + (1 + ^)" ~ a{a - l)(52/i-2, as \h\ -^ +oo. When a = 1 and 
1^1 > ^) li,jih.,6) reduces to 



7ij(/i, 5) = ^r7,,,5(/i) with B{h) =i{h~6) log ( 1 - - j + (/i + 5) log ( 1 + - 

Using the expansion of log(l ± x) as x — > leads to B{h) ~ 5^|/i|~^ as \h\ — )■ +oo, which 
implies the result. D 

Proposition H] and ([6]) lead to the following important remarks on the dependence struc- 
ture. For i y^ j and Hi + Hj ^ 1 : 

• If the two fractional Gaussian noises are short-range dependent (i.e. H^ < 1/2 and 
Hj < 1/2) then they are either short-range interdependent if pij 7^ or rjij 7^ 0, or 
independent if pij = rji^j = 0. 

• If the two fractional Gaussian noises are long-range dependent (i.e. Hi > 1/2 and 
Hj > 1/2) then they are either long-range interdependent if pij 7^ or rjij ^ 0, or 
independent if pij = Tjij = 0. This confirms the dichotomy principle observed in 

m- 

• In the other cases, the two fractional Gaussian noises can be short-range interde- 
pendent if Pij 7^ or rjij 7^ and Hi + Hj < 1, long-range interdependent if Pij 7^ 
or r]ij 7^ and Hi + Hj > 1 or independent if pij = rjij = 0. 

Moreover, note that when Hi + Hj = 1, whatever the nature of the two fractional Gaussian 
noises (i.e. short-range or long-range dependent, or even independent), they are either 
long-range interdependent if fjij 7^ or independent if fjij = 0. 

The following result characterizes the spectral nature of the increments of a mfBm. 
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Proposition 5 (Coeurjolly et al. [5]). Let Sij{-,6) be the (cross)-spectral density of the 
increments of size 6 of the components i and j , i.e. the Fourier transform of'-fij{-,6) 

Sijiu, 6) = ^ f e-'^-^,,{h, 5) dh =: Fr(7,,(-, 6)). 

(z) For all i,j and for all Hi, Hj, we have 

(9) S,,{u;, 6) = "-^ViH, + H, + 1)^^^^ x r,,(sign(c.)), 

TT \UJ\ ^' 3^ 

where 
(10) 

r (sisn(u)) = I P'^' '''' (^^^* + ^'^^ ~ i^Msign(a;) cos (f (if, + H,)) if H, + H, ^ 1, 
''^y ^ y >> \ p.^^ - i^fj^^jsign{uj) ifH, + Hj = 1. 

{ii) For any fixed 6, when Hi + Hj ^ 1 then we have, as w — )■ 0, 

(11) 

\Si,j{^,0)\r^^;^T{H, + Hj + l)d \u;\H.+H,-l " 

(iii) Moreover, when Hi + Hj ^ 1, the coherence function between the two components i 
and j satisfies, for all u 



(12) 



Si^i{u,6)Sjj{uj,6) 

T{H, + H, + 1)2 pI sin (f (if, + H,)Y + ^^^ cos (f (ii, + H,))' 



T{2Hi + l)T{2Hj + 1) sin(7riii) sin(7riij) 



(iv) When H^+Hj = 1, (E]) and (^ hold, replacing plj sin (f (ii^ + Hj)Y+7]lj cos (f (ii^ + Hj))' 

hPlj + TVlj- 

Proof. Tlie proof is essentially based on the fact that in the generalized function sense, for 
a > —1, 



-o-l 



Fr(|/i|") = --r{a + 1) sin (^a]\uj\ 

FT{hl) = — r(a + l)e-''^s"(^)5(°+^)|a;r°-\ 

FT(/i^) = — r(a + l)e''''snMt("+^)|a;|'"-\ 
27r 

FTihlog\h\) = i-^^^"^") 



2^2 
See |5] for more details. D 



6 AMBLARD, COEURJOLLY, LAVANCIER, AND PHILIPPE 

Remark 4. From this proposition, we retrieve the same properties of dependence and 
interdependence of Xj and Xj as stated after Proposition ^ 

2.3. Time reversibility. A stochastic process is said to be time reversible if X{t) = 
X{—t) for all t. As shows in |12], this is equivalent for zero-mean multivariate Gaussian 
stationary processes to KXi{t)Xj{s) = KXi{s)Xj{t) for s, t € M or that the cross covariance 
of the increments satisfies jijlhjS) = '~fij{—h,6) for h E M.. The following proposition 
characterizes this property. 

Proposition 6. A mfBm is time reversible if and only if r]ij = (or fjij = 0) for all 

Proof. If rji^j = (or fjij = 0), 'jij{h,6) is proportional to the covariance of a fractional 
Gaussian noise with Hurst parameter {Hi + Hj)/2 and is therefore symmetric. Let us prove 



the 


converse, 


. Let 


a - 


= Hi + Hj, then 


li 


.,(h,S)- 


7i,i(- 


-h, 


, 6) = aiffj X 






1 


-Vi,i 


^s: 


ign(/i — S) h - 


-6 



2sign(/i)|/i|° - sign(/i + 5)\h + 6\") if a ^ 1, 
\ fjij {{h - 6) log \h-6\- 2/ilog \h\ + {h + 6) log \h + 6\) if a = 1. 

Assuming 'Jijih, 6) — '-fij{—h, 6) equals zero for all h leads to rjij = (or r/j j = 0). D 

Remark 5. This result can also be viewed from a spectral point view. The time reversibility 
of a mfBm is equivalent to the fact that the spectral density matrix is real. Using (^, this 
implies ijij = (or fjij = 0). 

3. Integral representation 

3.1. Spectral representation. The following proposition contains the spectral represen- 
tation of mfBm. This representation will be especially useful to obtain a condition easy to 
verify which ensures that the functions defined by ([3]) and (j3]) are covariance functions. 

Theorem 7 (Didier and Pipiras, [I2]). Let {X(t))te_R be a mfBm with parameter {Hi, ■ ■ ■ , Hp) 
(0, 1)^. Then there exists a p x p complex matrix A such that each component admits the 
following representation 

P /• \tx "1 

(13) X,{t) = J2 / ^^^(A,,x;^'+V2 + A,,x-_''''-'^')B,{ dx), 

where for all j = l,...,p, Bj is a Gaussian complex measure such that Bj = Bj^i + 
iBj2 with Bji{x) = Bj^i{—x), Bj^2{x) = —Bj2{x), Bji and Bj2 are independent and 
E{Bj^i{dx)Bj^i{dxy) = dx, 1 = 1,2. 

Conversely, any p -multivariate process satisfying (IT3|) is a mfBm process. 



Proof. This representation is deduced from the general spectral representation of operator 
fractional Brownian motions obtained in [12]. By denoting — EI + 1/2 := diag(— iiTi + 
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1/2, ■ • ■ , —Hp + 1/2) we have indeed 



(14) X{t) = / '-^{x-^^"'A + x-''^"'A)B{ dx). 



D 



Any mfBm having representation (1131) has a covariance function as in Proposition [31 
The coefficients p, j, r^j j, pij and fjij involved in ([2]) and (jl]) satisfy 

(15) {AA*\, = ^r{H, + Hj + l)r,,(l), 

where Tij is given in (fTOj) and where A* is the transpose matrix of A. This relation is 
obtained by identification of the spectral matrix of the increments deduced on the one 
hand from ( IT3l) and provided on the other hand in Proposition [5l 

Given flT3l) . relation ( ITSl) provides easily the coefficients pij, rjij, pij and fjij which define 
the covariance function. The converse is more difficult to obtain. Given a covariance 
function as in Proposition |3l obtaining the explicit representation flT3|) requires finding a 
matrix A such that (fT5|) holds. This choice is possible if and only if the matrix on the right 
hand side of f flSj) is positive semidefinite. Then a matrix A (which is not unique) may be 
deduced by the Cholesky decomposition. When p = 2, an explicit solution is the matrix 
with entries, for i,j = 1, 2, 



^i,j ^i,i 



pij sin {-{Hi + Hj)) + r]ijj "'-^ cos y-^i^i 





sin (^^{H, + Hj)j - r]ij cos (^^{Hi 




where \ij = — ^-= — '' and Cj ,■ is given in (fT2l) . provided Hi + H2 ^ 1. 

2V7r ^r(2ifj + 1) sin{Hj7r) 

When H1 + H2 = 1, the same solution holds, replacing pij by pij and r]ij cos (f (-f^i + -^^2)) 
by -lVi,r 

3.2. Moving average representation. In the next proposition, we give an alternative 
characterization of the mfBm from an integral representation in the time domain (or moving 
average representation). 

Theorem 8 (Didier and Pipiras, [I2]). Let {X{t))t(zR be a mfBm with param,eter {Hi, ■ ■ ■ , Hp) G 
(0, 1)^. Assume that for alii G {l,...,p}, Hi 7^ 1/2. Then there exist M~^,M~ two p x p 
real matrices such that each component admits the following representation 
(16) 



X^it) = X^ / M+ ((t - x)^^-' - (-x)?-^) + M,- ((t - x)^- 
,=1 J^ 



x)_^-^) Wj{dx), 
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with W{dx) = {Wi{dx), ■ ■ ■ ,Wp{dx)) is a Gaussian white noise with zero mean, indepen- 
dent components and covariance E,Wi{dx)Wj{dx) = Sijdx. 

Conversely, any p-multivariate process satisfying ( IT6|) is a mfBm process. 

Proof. This representation is deduced from the general representation obtained in [12]. D 

Remark 6. When Hi = 1/2 for each i G {1, ...,p}, it is shown in jl2] that each component 
of the mfBm admits the following representation : 

P r. 

Xi{t) = XI / ^i^j(sign(t - x) - sign(x)) + M". (log \t - x\ - log |x|) Wj{dx). 

Our conjecture is that this representation remains valid when Hi = 1/2 whatever the values 
of other parameters Hj, j ^ i. 

Starting from the moving average representation flT6|) . using results in j2^, we can specify 
the coefficients pij, rjij, pij and fjij involved in the covariances ([3]) and (jll) (see 120]). More 
precisely, let us denote 

M+(M+)' =(«++) , M-iM-y={arr), M+(M-)' = («+r) 

where M' is the transpose matrix of M. The variance of each component is equal to 

2 B{Hi + .5,Hi + .5) r ++ , — ,„ X +_^ 

where B{-, ■) denotes the Beta function. 
Moreover, if Hi + Hj ^ 1 then 

B{Hi + .5, Hj + .5 ) 

sm{{Hi + Hj)n) 
{«/ + «^7)(cos(if.7^) + cos(if,7r)) - (a+r + «"+) sm{{H, + H^)7r)} 
B{H, + .5, Hj + .5 ) 

sm{{H, + H,)n) 

{«/ - «.7)(cos(i:r,7r) - cos(i:r,7r)) - («+" - «"+) sin((i:r, + H,)7r)} . 

If iJi + Hj = l then 

„,^^ ^ ^^ ^, rsin(i^j7r) + sin(iJo7r) , ,, , , ,^ 

a,a,p,, = B{H, + .5, H^ + .5)| ^ ' ' ^ ^^^K/ + «."/) " <" " " j|' 

Conversely, given a covariance function as in Proposition |3l if Hi ^ 1/2 for all i, one 
may find matrices M"*" and M~ such that (IT6|) holds, provided the matrix on the right 
hand side of ( TT5|) is positive semidefinite. Indeed, in this case, a matrix A which solves 
(IT5|) may be found by the Cholesky decomposition, then M"*" and M~ are deduced from 
relation (3.20) in [I2]: 



O-iO-jpiJ = _,_,frT . TT\_\ ^ 



(^iO-jVid = _.___UTT . /^^_^ ^ 
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M± = y|(Dr%±/^2"%) 



where A = Ai + 1^2 and 



Di = diag sm{ivHi)TiHi + -),..., sm(7ri/p)r(i7p 



-),..., sm{7iHp)T{Hp+- 



D2 = diag (cos{7rHi)T{Hi + ^), • • • , cosiTTHp)T{Hp + ^) j . 

3.3. Two particular examples. Let us focus on two particular examples which are quite 
natural: the causal mfBm (M~ = 0) and the well-balanced mfBm (M~ = M+). In the 
causal case, the integral representation is a direct generalization of the integral representa- 
tion of Mandelbrot and Van Ness j22] to the multivariate case. The well-balanced case is 
studied by Stoev and Taqqu in one dimension [27]. With the notation of the two previous 
sections, we note that the causal case (resp. well-balanced case) leads to Ai = tan{7rH)A2 
(resp. A2 = 0), where tan(7riJ) := diag (tan(7riJi), . . . , tan(7rifp)). In these two cases, 
the covariance only depends on one parameter, for instance pij (or Pij). Indeed we easily 
deduce rjij (or fjij) as follows : 

• in the causal case i.e. M~ = or equivalently Ai = tan{'n'H)A2 : 

Vi,j = -Pi,j tan(^(i/i + Hj)) tan{^{Hi - Hj)) if H, + H^ ^ 1, 

• in the well-balanced case i.e. M~ = M^ or equivalently ^2 = : 

7],j = iiH, + H,^l, 
fjij = iiHi + Hj = 1. 



Remark 7. From Proposition\^ the well-balanced mfBm is therefore time reversible. 

3.4. Existence of the covariance of the mfBm. In this paragraph, we highlight some 
of the previous results in order to exhibit the sets of the possible parameters [pij,riij) or 
{Pi,j^Vi,j) ensuring the existence of the covariance of the mfBm. 

For i,j = l,...,p, let us give {Hi,Hj) G (0,1)^, (crj,crj) G IR+ x M+ and (p^j rjij) e 
[-1, 1] X M with pj^i = pij and r]j^i = -r]ij if Hi + Hj ^ 1, or (p^j, ?7ij) G [-1, 1] x M with 
Pj.i = Pi,j and fjj^i = -fjij if Hi + Hj = 1. 

For this set of parameters, let us define the matrix S(s,t) = (Sjj(s,t)) as follows : 
Sj^j(s,t) is given by ([2]) and Sjj(s,t) is given by ([3]) when Hi + Hj ^ 1 and (|4]) when 
Hi + H, = 1. 
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Proposition 9. The matrix S(s,t) is a covariance matrix function if and only if the 
Hermitian matrix Q = {T{Hi + Hj + l)rjj(l)) with Tij defined in f ITOl) . is positive semidef- 
inite. When p = 2, this condition reduces to Ci 2 < 1 where Ci 2 is the coherence defined 

bym- 

Proof. First, note that since pji = pij and rjj^i = —rjij , Q is a Hermitian matrix. Now, if 
Q is positive semidefinite, then so is the matrix (27r)~^(o"j(Tj(5ij). Therefore there exists a 
matrix A satisfying (fTSj) . From Theorem[71 there exists a mfBm having S(s, t) as covariance 
matrix function. Conversely, if T,{s,t) is a covariance matrix function of a mfBm then the 
representation (IT3ll holds and by ( IT5|) . the matrix Q is positive semidefinite. 

When p = 2, the result comes from the fact that Q is positive semidefinite if and only if 
det{Q) > or equivalently Ci^2 < 1- CH 



When p = 2, for fixed values of {Hi,H2) the condition Ci^2 < 1 means that the set of 
possible parameters (pi,2,''7i,2) is the interior of an ellipse. These sets are represented in 
Figure [1] according to different values of Hi and H2. Note that, in order to compare the 
cases H1+H2 7^ 1 and H1+H2 = 1, we have reparameterized 771,2 by ri[ 2 := f]i,2{^—Hi—H2). 
In such a way, the second ellipse becomes the limit of the first one as Hi + H2 ^ 1 (see 
also Remark [2]). 

Let us underline that the maximum possible correlation between two fBm's is obtained 
when ?7i 2 = 0, i.e. when the 2-mfBm is time reversible according to Proposition O 

Remark 8. When Hi = . . . = Hp = H ^ 1/2, the matrix Q rewrites Qij = T{2H + 
l){sin{'iTH)pij — ir]ij cos{7tH)) and 

• if the mfBm is time reversible, i.e. rjij = (for i,j = l,...,p), then Q is a 
correlation matrix and is therefore positive- semidefinite for any \pij\ < 1, 

• when p = 2, the set of possible values for {pi^2,Vi,2) associated to H and 1 — H are 
the same. 

In the particular case of the causal or the well-balanced mfBm, the matrix S(s,t) can 
be expressed through the sole parameter pij. Let X{Hi,H2) the function which equals to 
cos(|(iJi — H2)Y ii^ the causal case and which equals 1 in the well-balanced case. The 
maximal possible correlation when p = 2 is given by 

_ V{2Hi + l)V{2H2 + 1) sin(7ri7i) sin(7rff2) .. , , ^ ^. 
^^'^ " T{Hi + H2 + lf sHl{Hi + H2)f "" ^^^1'^^^- 

Figure [2] represents |pi.2| with respect to {Hi,H2). 

Figures [1] and [2] illustrate the main limitation of the mfBm model. Under self-similarity 
condition ([1]), it is not possible to construct arbitrary correlated fractional Brownian mo- 
tions. For example, when Hi = 0.1 and H2 = 0.8, the correlation cannot exceed 0.514. 
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HI =0.1, H2=0.1 
HI =0.1, H2=0.2 
HI =0.1, H2=0.3 
HI =0.1, H2=0.4 
HI =0.1, H2=0.5 



HI =0.1, H2=0.6 
HI =0.1, H2=0.7 
HI =0.1, H2=0.£ 
HI =0.1, H2=0.9 






Figure 1. Various examples of possible values for {pi,2,Vi2) with 77^3 •= 
771,2(1 - Hi - H2) when Hi + H2 j^ I and (pi,2,''?i,2) when Hi + H2 = 1, 
ensuring that S(s,t) is a covariance matrix function in the particular case 
p = 2. 



4. The mfBm as a limiting process. 

A natural way to construct self-similar processes is through limits of stochastic processes. 
In dimension one, the result is due to Lamperti [18]. In [16], an extension to operator self- 
similar processes is given. A similar result for the mfBm is deduced and stated below. In 
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Maximal values of p - Causal case 



Maximal values of p - Well-balanced case 





Figure 2. Maximal values of the absolute possible correlation parameter |/9i,2| 
ensuring that S(s,i) is a covariance matrix function in the case p = 2, in terms of 
Hi and H2 for the causal and well-balanced mfBm. 

the following, a ^-multivariate process (X(t))tgiK is said proper if, for each t, the law of 
X{t) is not contained in a proper subspace of M^. 

Theorem 10. Let (X(t))tgR be a p-multivariate proper process, continuous in probability. 
If there exist a p-multivariate process (Y(t))t(zR and p real functions ai,....,ap such that 

(17) (aiHFi(nt), . . . , ap{n)Y,{nt)) ^^^^ X{t), 

fidi 

then the multivariate process {X(t)) is self-similar. Conversely, any multivariate self- 
similar process can be obtained as a such limit. 



Proof. The proof is similar to Theorem 5 in (TB]. Fix A; G N and r > 0. For each T G M'^ we 
denote X(T) := (X(Ti), . . . , X(Tk)). Let Vr^k be the set of all invertible diagonal matrices 
a such that, for all T eR'', X{rT) = aX{T). 

Let us first show that P^.fc is not empty. According to flT7|) . we have 



and 



diag(ai(n), . . . , ap{n))Y{nrT) > X{rT), 



diag(ai(rn), . . . , ap{rn))Y (nrT) > ^(T). 

d 

Since {X{t)) is proper, diag(ai(n), . . . , ap{n)) and diag(ai (rn), . . . , ap{rn)) are invertible 
for n large enough. Then, Theorem 2.3 in [21] ensures that a„ defined by 

a„ = diag(ai(?7,), . . . , ap(ra))diag(ai(nr), . . . , ap{nr))~^ 

has a limit in Pr-,fc- Moreover if a is a limit of a„ then X{rT) = aX(T) and thus P^.fc 7^ 0- 
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It is then straightforward to adapt Lemma 7.2-7.5 in |16| for the subgroup Pr,fc; which 
yields that for each r, flfcPr.fc is not empty. Therefore, for any fixed r > 0, there ex- 
ists a G (Ik'Dr^k such that {X{rt)) and {aX(T)) have the same finite dimensional dis- 
tributions. Theorem 1 in [16] ensures that there exists {Hi, . . . ,Hp) G (0,1)^ such that 
a = diag(r^i, . . . , r^p). The converse is trivial. 

D 

As an illustration of Theorem [101 the mfBni can be obtained as the weak limit of partial 
sums of sum of linear processes (also called superlinear processes, see |33)). Some examples 
may be found in [8] and [1^. In Proposition [TT] below, we give a general convergence result 
which allows to reach almost any mfBm from such partial sums. The unique restriction 
concerns the particular case when at least one of the Hurst parameters is equal to 1/2. 

Let {ej{k))k&z, j = 1, ■ ■ ■ ,p he p independent i.i.d. sequences with zero mean and unit 
variance. Let us consider the superlinear processes 

p 

(18) z,{t) = J2Y1 ^m(^ - ^)^^(^)' ^ = 1. • • • .p. 

j=i fcez 

where ipij^k) are real coefficients with Ylikez'^l,]^^) ^ °*^- 

Moreover, we assume that ipij{k) = 'ip^jik) + 'ipijik) where ip^j^k) satisfies one of the 
following conditions: 

(i) ip^j{k) = [afj + o(l)) k^-'~ as |A;| -> oo, with < dfj < ^ and «+• ^ 0, 

d+ -i 
(ii) i^tj(^) = «i + o(l)) k^'' as \k\ -> oo, with -| < rf+. < 0, EfcgzV^^^'i(^) = ^nd 

(iii) EfcGZ |V^*ti(^)| < °° ^^d l^t «i •= Efcez^»\'j(^) ^ 0' «?+• := 0. 
Similarly, ip^jik) is assumed to satisfy (i), {ii) or {iii) where k^, dfj and afj are replaced 
by /;;_, d^j and a^j. 

Proposition 11. Let di = max{df^, d~^, ■ ■ ■ , dfp, d~p), for i = 1, . . . ,p. Consider the vector 
of partial sums, for t eW, 

([tit] [nr] 

n^^^~W2)J2z,{t),...,n-'p-^'/'^J2^p(^) 
t=i t=i 

Then the finite dimensional distributions of {Sn{T))reR converge in law towards a p-mfBm 
{X{r))rm- 

• When di ^ 0, {Xi{T))reR is defined through the integral representation [T^) where 
M^j = Oit,jd~^^d+^=d, (^nd M-j = a-.d-^l^-_^^. 

• When di = 0, Xi{r) = jyi=ii^i'Ad+ =o + ^^Ad" =o)^j{'^)> where Wj is a standard 
Brownian motion. 
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Moreover, if for all j = l,...,p, E{ej{0)^°') < oo with a > 1 V (1 + 2dmax)~^ where 
dmax = maxi{di}, then S'.„(.) converges towards the p-mfBm X{.) in the Skorohod space 
P([0,1]). 

Sketch of proof . We focus on the convergence in law of 5'„(r) to X(t), for a fixed r in M, 
the finite dimensional convergence is deduced in the same way. We set for simplicity r = 1. 
According to the Cramer- Wold device, for any vector (Ai, . . . , Xp) G MP, we must show 
that A'S'„(1) converges in law to A'X(l). We may rewrite A'5'„(l) as a sum of discrete 
stochastic integrals (see [29] and |1]) : 

p n 

i=l t=l 

p p n 

= E ^'^ E E ^""'"^'^'^ E(^i(^ -k)+ ^IM - k))eAk) 
4=1 j=i kez t=i 

(19) =E^''E /(/i» + /M»)^.>nW' 

where the stochastic measures Wj^n, j = ^, ■ ■ ■ iP are defined on finite intervals C by 

WaC)=n-'/' J2 e,{k), 

k/neC 

and where /^"t- „, /j~- „ are piecewise constant functions defined as follows: denoting [x"| the 
smallest integer not less than x, we have for all x G M, 



/i» = ^-''E<.(^-M)' 



i=l 



respectively /^ ,>(a;) = n '^^ Er=i^i,,(^ - M). 



The following lemma states the convergence of a linear combination of discrete stochastic 
integrals as in (fT9|) . A function is said n-simple if it takes a finite number a nonzero constant 
values on intervals {k/n, {k + 1)/?^], A; G Z. 

Lemma 12. Let (/i,n, ■ ■ ■ , fp,n)n&N be a sequence of p n-simple functions in L^(]R). // 
for any j = 1, . . . ,p, there exists fj G L^(M) such that J^^ \fj,n{x) — fj{x)\^dx — > 0, then 
Yl^=i lM.fj,nix)^j,n{dx) converges in law to ^^j=i J^ fj{x)Wj{dx) , where the Wj's are in- 
dependent standard Gaussian random measures. 

When p = 1, this lemma is proved in |28| . The case p = 2 is considered in [1] and the 
extension to p > 3 is straightforward. 

From Lemma [12] and flT9]) , it remains to show that 

2 



lim 

n— >oo 



a^ 



/S»-^((i-^)±-(-^)±)i.±. 



di ^'' '~'^ ' ~"'^' '"1,='^^ 



dx = 0, 
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where we agree that c?j^^((l — a^)± — (— a^)±) = l[o,i](2;) when di = 0. Below, we only consider 
the pointwise convergence of /j^-„(x), for x G M, when df, = di. The convergence in L^ is 
then deduced from the dominated convergence theorem (see [2S], [22], |1] for details). It 
also follows easily that, when d^, < di, f^ |/j •„(x)p(ia; — )■ 0. 

Under assumption (i), note that since di > 0, (1 — x)^!: — {~x)^ = di J {t — x)j!:~ dt. 
We have, for any x G M, 

n 
t=l 

= n-"^ J\^j{\t-] - \nx-\)dt 



= n-'^ / {at,+oil))i\t]-\nx])r'dt 
Jo 

= /' («£• + -(1)) (^^^^^i^)'"' dt -^ at, j\t - xrr'dt 

Under assumption (ii), di < 0. When a; < 0, (1 — x)^' — {—x)^ = di J^ {t — xY^~^dt and 
the convergence of f^j^ix) can be proved as above. When x > 1, (1 — x)^' — (— x)^' = 
= /+-„(x). When < x < 1, (1 - x)+' - (-x)J = -rfi/+°°(t - xY^'^dt and, since 
EfcezV^*ti(^) = 0' we have 

n n—\nx] 

Therefore, 



fi,ni^) = -n-'^ «. + o(l))([tl)'^-Vt 

>/n— [nx] 



~00 /" + 00 

-a+ / f'-'dt = -a+ / (t - xY'-'dt 

'l~x ' Jl 



This proves ftj^ni^) ~^ '^i "'^'^ij((l ~ ^)+ ~ ('~^) + )) ^^ ^'^Y a^ ^ IR, under assumption 
The same scheme may be used to prove that /j~„(x) — ?■ d'^^a'^Ml — x)J — (— x)j) under 
assumption (ii), noting that 

{0 when x < 0, 

-di J'^^it - xY'-^dt when < x < 1, 
di J^it - xf2~^dt when x > 1. 
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Under assumption (iii), 

n n— \nx] 

i=l t=l-lnx~\ 

Since ^^^^ i'^j{t) < oo, ftjni^) ^ for all x ^ [0, 1]. When x G [0, 1], we have ftjni^) ^ 

Therefore, the first claim of the theorem is proved, i.e. the convergence in law of the finite 
dimensional distribution of (5'„(T)),-giK to (X(r))T-gR. To extend this convergence to a func- 
tional convergence in 'D{[0, 1]), it remains to show tightness of the sequence (S'„(r))T-g[o,i]. 
This follows exactly from the same arguments as in the proof of Theorem 1.2 in |3|. D 

5. Synthesis of the mfBm 

5.1. Introduction. The exact simulation of the fractional Brownian motion has been a 
question of great interest in the nineties. This may be done by generating a sample path of 
a fractional Gaussian noise. An important step towards efficient simulation was obtained 
after the work of Wood and Chan |32] about the simulation of arbitrary stationary Gaussian 
sequences with prescribed covariance function. The technique relies upon the embedding 
of the covariance matrix into a circulant matrix, a square root of which is easily calculated 
using the discrete Fourier transform. This leads to a very efficient algorithm, both in terms 
of computation time and storage needs. Wood and Chan method is an exact simulation 
method provided that the circulant matrix is semidefinite positive, a property that is not 
always satisfied. However, for the fractional Gaussian noise, it can be proved that the 
circulant matrix is definite positive for all H G (0, 1), see j9l [13] . 

In |7], Wood and Chan extended their method and provided a more general algorithm 
adapted to multivariate stationary Gaussian processes. The main characteristic of this 
method is that if a certain condition for a familiy of Hermitian matrices holds then the 
algorithm is exact in principle, i.e. the simulated data have the true covariance. We present 
hereafter the main ideas, briefiy describe the algorithm and propose some examples. 

Remark 9. Other approaches could have been undertaken (see |3] for a review in the 
case p = 1). Approximate simulations can be done by discretizing the moving- average or 
spectral stochastic integrals l[TS\) or D^). [B] also proposed an approximate method based on 
the spectral density matrix of the increments for synthesizing multivariate Gaussian time 
series. Thanks to Proposition\^ this could also be envisaged for the mfBm. 

5.2. Method and algorithm. For two arbitrary matrices A = (Aj,fc) and B, we use A^B 
to denote the Kronecker product of A and B that is the block matrix (AjkB). 

Let AX := AiX denotes the increments of size 1 {6 = 1) of a mfBm. We have 
AX = (AX(t))4gK = ((AXi(t), . . . , AXp(t))')te]R. The aim is to simulate a realization 
of a multivariate fractional Gaussian noise discretized at times j = 1, . . . , n, that is a real- 
ization of (AX(1), . . . , AX(n)). Then a realization of the discretized mfBm will be easily 
obtained. 
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We denote by AX^") the merged vector AX^") = (AX(1)', . . . , AX(n)')' and by G its 
covariance matrix. G is the npxnp Toephtz block matrix G = (Gdi— j|))jj=i^...^„ where for 
h = 0, . . . , n — 1, G{h) is the pxp matrix given by G{h) := {ij^kih))- ^^-^ . The simulation 
problem can be viewed as the generation of a random vector following a A/'„p(0, G). This 
may be done by computing G^^^ but the complexity of such a procedure is 0{p^n^) for 
block Toeplitz matrices. To overcome this numerical cost, the idea is to embed G into the 
block circulant matrix C = circ{C(j), j = 0, . . . , -m — 1}, where m is a power of 2 greater 
than 2(n — 1) and where each C{j) is the pxp matrix defined by 

f G{3) if < J < m/2 

(20) C{j)=l i(G(j) + G(j)') if J = m/2 

I G(j — m) if m/2 < j < m — 1. 

Such a definition ensures that C is a symmetric matrix with nested block circulant structure 
and that G = {C{j),j = 0, ... ,n — 1} is a submatrix of C. Therefore, the simulation of 
a A/'„p(0,G) may be achieved by taking the n "first" components of a vector Mmp{0,C), 
which is done by computing C^^^. The last problem is more simple since one may exploit 
the circulant characteristic of G: there exist m Hermitian matrices B{j) of size pxp such 
that the following decomposition holds 

(21) G = {J® Q diag(5(j), J = 0, . . . , m - 1) ( J* ® Jp), 

where Q is the m x m unitary matrix defined for j. A; = 0, m — 1 by Jj^ = e"^''^-''^'/™. The 
computation of C^'^ is much less expensive than the computation of G^^^ since, as in the 
one-dimensional case {p = 1), (12T]) will allow us to make use of the Fast Fourier Transform 
(FFT) which considerable reduces the complexity. 

Now, the algorithm proposed by Wood and Chan may be described through the following 
steps. Let ?72 be a power of 2 greater than 2{n — 1). 
Step 1. For 1 < u < v < p calculate for /c = 0, . . . , ?Ti — 1 

m—l 

5«,.(A;) = 5^a,.(j>-''"^'/'" 
where Gu,v{j) if the element {u,v) of the matrix G{j) defined by (120|) and set By^u{k) = 

Step 2. For each j = 0, . . . , m — 1 determine a unitary matrix R{j) and real numbers S,u{j) 
(w = l, • • • ,P) such that B{j) = R{j) diag(ei(j), • • • , Uj)) RUY- 

Step 3. Assume that the eigenvalues ^i(j), . . . ,^p(j) are non-negative (see Remark [TT]) and 
define 5(j) = RU) dmgi^/^j), ..., ./^j)) R{jy. 
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Step 4. For j = 0, . . . , m/2 generate independent vectors f/(j), V{j) ~ A/^(0, /) and define 

1 ^ / v^f/(j) for J = 0, 



^^^ v/2^ I f/(j) + i^(j) for J = 1, . . . , f - 1, 
let Z{m - j) = Z{j) for j = f + 1, . . . , m - 1 and set W{j) := B{j)Z{j). 

Step 5. For u = 1, . . . ,p calculate for /c = 0, . . . , rra — 1 

m—l 

and return |AX„(A;), l<n<j», A; = 0,...,n — l}. 

Step 6. For u = 1, . . . ,p take the cumulative sums AXu to get the u — th component X„ 
of a sample path of a nifBni. 

Figure [3] gives some examples of sample paths of mfBm's simulated with this algorithm. 

Remark 10. Let us discuss the computation cost of the most expensive steps, that is steps 
1, 2 and 5. Step 1 requires Hi^^t_i applications of the FFT of signals of length m, Step 2 
needs m diagonalisations ofpxp Hermitian matrices and Step 5 requires p applications of 
the FFT of signals of length m. Therefore, the total cost, K{m,p) equals 

/ p(p + 1) \ 

K{m,p) = O I TTilogm J + 0{mp^) + 0{pm\ogm). 



Remark 11. The crucial point of the previous algorithm lies in the non-negativity of the 
eigenvalues C,i{j), • • • , ^p{j) for any j = 0, . . . ,m — 1. In the one-dimensional case (when 
p = 1) Steps 2 and 3 disappear, and in Step 1, Buik) corresponds to the k—th eigenvalue 
of the circulant matrix Cu with first line defined by Cii(j) = 7ii(j) for < j < m./2 and 
'~f{m — j) for j = m/2 + 1, . . . , ??7. — 1. For the fractional Gaussian noise, it has been proved 
by Craigmile [H] for H < 1/2, and by Dietrich and Newsam jT3] for H > 1/2 that such 
a matrix is semidefinite-positive for any m (and so for the first power of 2 greater than 
2{n — l)). In the more general casep > 1, the problem is much more complex: the quantities 
Bu,v{k) are not necessarily real, and the establishment of a condition of positivity for the 
matrix Bu,v{k) does not seem obvious. When the condition in Step 3 does not hold. Wood 
and Chan suggest to either increase the value of m and restart Steps 1,2 or to truncate the 
negative eigenvalues to zero which leads to an approximate procedure. These problems are 
not addressed in this paper. Let us assert that for the simulation examples presented in 
Figure\^ we have observed that this condition is satisfied for m equal to the first power of 
2. 
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Figure 3. Examples of sample paths of length n = 1000 normalized causal 
(top), well-balanced (middle) and general (bottom) multivariate fractional 
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for these different choices of parameters. For convenience, the sample paths 
of the left column have been artificially shifted for visibility. 
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